
#include "functions.h"
#include <Eigen/Dense>
#include <cassert>
#include <iostream>
#include <vector>
using namespace std;
using namespace Eigen;

int main() {

  constexpr double re_num{1000000.0};
  VarArray<Point, 2, 0> grid;
  string grid_file_name{"flate_test.x"};
  ReadMeshASSCII(grid_file_name, grid);

  int npoint_in_ydir = grid.Dim(1);
  int npoint_in_xdir = grid.Dim(0);
  //   cout << npoint_in_ydir << endl;
  // 注意
  VarArray<Vector2d, 2, 0> flow(npoint_in_xdir, npoint_in_ydir);
  InitializeFlow(flow);

  VarArray<Matrix2d, 2, 0> jacob_matrix(3, npoint_in_ydir);
  vector<Vector2d> fx_k(npoint_in_ydir);
  vector<Vector2d> residual(npoint_in_ydir);
  //  迭代求解
  for (int i = 1; i < npoint_in_xdir; ++i) {

    int count = 0;
    double tmp{0.0}, res1, resn;
    do {

      count++;
      // // 计算之前未知时间层需要使用已知时间层初始化
      InitializeNextLine(flow, i);
      CalcJacobianMatrix(flow, i, grid, jacob_matrix, re_num);

      CalcFXk(flow, i, grid, re_num, fx_k);

      CalcResidual(jacob_matrix, fx_k, residual, grid);

      UpdateFlow(flow, residual, i);

      resn = CalcNormal(residual);
      if (count == 1)
        res1 = resn + 2.0e-24;
      tmp = log10(resn / res1);
    } while (tmp > -5.0 && resn > 1.0e-6 && count < 80);

    cout << count << endl;
    cout << "i = " << i << endl;
    exit(EXIT_FAILURE);
  }

  cout << jacob_matrix(0, 1) << endl;
  return 0;
}